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We present a model of the hydraulic fracturing of heterogeneous poroelastic media. The formalism 
is an effective continuum model that captures the coupled dynamics of the fluid pressure and the 
fractured rock matrix and models both the tensile and shear failure of the rock. As an application 
of the formalism, we study the geomechanical stress interaction between two injection points during 
hydraulic fracturing (hydrofracking) and how this interaction influences the fracturing process. For 
injection points that are separated by less than a critical correlation length, we find that the frac- 
turing process around each point is strongly correlated with the position of the neighboring point. 
The magnitude of the correlation length depends on the degree of heterogeneity of the rock and is 
on the order of 30 — 45 m for rocks with low permeabilities. In the strongly correlated regime, we 
predict a novel effective fracture-force that attracts the fractures toward the neighboring injection 
point. 



I. INTRODUCTION 

Hydrofracking is a technology that utilizes highly pres- 
surized fluid to create fracture networks in rock layers 
with low permeabilities. A fracking fluid is injected into 
a cased wellbore, and the parts of the reservoir to be 
fractured are accessed by perforating the casing at the 
correct locations (Fig. ^p)- The injection well can be 
vertical or horizontal, and several injection points may 
be active during the hydraulic stimulation of the system 
(Fig. [l^-b). The highly pressurized fluid that flows into 
the reservoir increases the geomechanical stress around 
the injection points and causes the rock to fracture. Ide- 
ally, the hydrofracking creates long, distributed cracks 
that connect a large area of the reservoir to the well. A 
comprehensive understanding of this fracturing process 
is therefore crucial to optimize the functionality of the 
reservoir. 

The creation of a fracture changes the geomechanical 
strain energy of the system. If the reservoir contains sev- 
eral fractures, the modification of the strain energy me- 
diates an effective interaction between the cracks. The 
interaction of fractures has been studied theoretically in 
several works using bo th dir ect numerical simulations^^ 
and statistical method^^^^. A geomechanical stress in- 
teraction also exists between two injection points during 
hydrofracking. A thorough investigation of this interac- 
tion and its consequences for the fracturing process is 
lacking. 

The failure of homogeneous materials is well under- 
stood. A fracture is created when the strain energy re- 
lease rate exceeds a critical value, and it propagates along 
the direction that maximizes the energy releas^^. The 
failure of heterogeneous materials, in contrast, depends 
on many microscopic mechanisms and is challenging to 
modeP. To understand the fracturing process of a het- 
erogeneous material on the micro-scale, statistical mod- 
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FIG. 1: (a) Hydrofracking using two vertical wells, (b) Hy- 
drofracking using a perforated horizontal well, (c) A sim- 
plified, discretized 2D representation of the reservoirs shown 
in (a) and (b). (d) A heterogeneous system is modeled by 
distributing the material strengths of the discrete volume el- 
ements according to the Weibull distribution. The degree of 
disorder is tuned by varying the Weibull modulus, m. 



els such as the fuse - and the beam- models have been 
particularly usefuJ^. In these models, the heterogeneity 
of the system is taken into account by defining a local, 
position-dependent material strength that is drawn from 
a statistical distribution. The fracturing process is there- 
fore not entirely controlled by the energy release rate, but 
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also depends on the disorder of the system and the dis- 
tribution of the local strengths. Depending on the degree 
of disorder, the statistical models have shown that het- 
erogeneous systems exhibit fracturing regimes that are 
distinct from the single-crack development observed in 
homogeneous systemJ^. 

Several models for hydrofracking have been devel- 
opecffl, but a weakness of most of these models is that 
they do not account for the heterogeneities of the rock. 
However, some research based on extensions of the sta- 
tistical beam model has been applied to the study of 
hydrofracking of heterogeneous system^SEltlll^ Recently, 
these studies have been further developed to provide a 
better description of the rock matrix. In Wangen and 
Wang et a/.^, the dynamics of the rock matrix are mod- 
eled using poroelastic theory, and the heterogeneity is 
treated by distributing the local material strength ac- 
cording to a probability distribution. 

In the present paper, we develop a model of the hy- 
drofracking of heterogeneous, poroelastic media. The for- 
malism captures the coupled dynamics of the fractured 
rock matrix and the fluid pressure. The main parts of the 
model, which describes the fracturing event and the ef- 
fective coupling between the pressure and the fractures, 
are based on a previous study by Wang en '^^ . We ex- 
tend this work by taking into account the anisotropic 
fluid flow and the shear failure of the material. Hetero- 
geneities are treated statistically by distributing the local 
strength of the material according to the Weibull distri- 
bution (Fig.[l]i). We apply the formalism to the study of 
the mechanical stress interaction between two injection 
points during hydrofracking, and investigate the depen- 
dency of the interaction on the disorder of the rock. For 
two points that are separated by a distance that is smaller 
than a critical correlation length, we flnd that the fractur- 
ing process around each injection point is strongly corre- 
lated with the position of the neighboring point. The crit- 
ical correlation length at which this strongly correlated 
regime occurs depends on the degree of heterogeneity, 
with correlation lengths of approximately 20 m for highly 
disordered systems and 45 m for weakly disordered sys- 
tems. In the correlated regime, we predict a novel effec- 
tive fracture force that attracts the fracture toward the 
neighboring injection point. Our results are important 
for optimizing the hydraulic stimulation of reservoirs. For 
well perforations that are separated by a distance that is 
less than the critical correlation length, the results im- 
ply a reduced effect of the stimulation because the frac- 
tures are attracted toward neighboring injection points. 
Knowing the correlation length of the system is therefore 
crucial for creating an effective and long-ranging fracture 
network. 

This paper is organized in the following manner. In 
Sec. [Hj we present the theory and governing equations 
of our hydrofracking model. The section concludes with 



a pseudo code of the algorithm. Sec. |ni| describes the 
numerical solution strategy. The next two sections ap- 
ply the formalism to the study of the interaction of two 



injection points during hydrofracking. Sec. [TV| provides 
a description of the model system that we consider, and 
Sec. [V| presents our flndings. We conclude and summa- 
rize our results in Sec. IVII 



II. GOVERNING EQUATIONS 

In this section, we consider a poroelastic system and 
develop the mathematical description that captures the 
coupled dynamics of the fractured rock matrix and the 
fluid. At the end of the section, the tensile and shear 
failure criteria are deflned. 



A. Poroelastic Theory 

Under the action of applied forces, an elastic medium 
exhibits deformations. The deformations can be in the 
form of a change in the shape of the object (without a 
change in its volume), referred to as a shear deformation, 
or a compression or stretching, referred to as a volumet- 
ric deformation. Let r denote the position of a material 
point in the medium before the deformation, and let 
denote its value after the deformation. The displace- 
ment of the point is then given by the displacement fleld 
u(r) = — r. The displacement fleld is a function of 
the position r, i.e., it describes how each material point 
moves under a deformation. For small displacements, the 
state of the elastic system is completely described by the 
strain tensoi^^: 



dui . du 



'j_ 

dxj 



(1) 



The diagonal elements of the strain tensor represent local 
volumetric deformations around each point r, while the 
off-diagonal terms capture the shear deformations. In 
mathematical terms, the displacement fleld is a vector 
fleld, and the strain tensor is a tensor fleld over the space 
formed by the elastic object. 

The dynamics of the elastic medium are described by 
the following fundamental law of elasticity Xheoiy^^: 



' dt^ 



dxk ' 



(2) 



Here, p is the mass density, Fi represents the external 
forces, and cr^j is the stress tensor arising from the in- 
ternal stresses. The internal stresses are caused by the 
intermolecular forces that occur because of the relative 
displacement of the molecules under the deformation. In 
Eq. Q (and in what follows), we apply the Einstein sum- 
mation convention and sum over repeated indices. Usu- 
ally, the only external force that appears is the gravita- 
tional force, F = pg, where g is the gravitational acceler- 
ation. The equilibrium state of the system is determined 
by solving the stationary equation, i.e., Eq. ([2| with 
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d^Uijdt^ = 0. The elastic energy stored in a strained 
system is-^-^: 



E 



= J a. 



(3) 



In this study, we concentrate on a porous, elastic rock 
system with water-filled pores. In the linear response 
regime for an isotropic, poroelastic system, the stress ten- 
sor can be written phenomenologically aa^ 



where 



= 



2Ge, 



(4) 



(5) 

The tensor a[j is the Terzaghi effective stress tensor, 
which describes the stresses that act only on the rock ma- 
trix, while aij represents the stresses acting on the total 
fiuid-rock system; is the trace of the strain tensor, and 
A > and G > are elasticity coefficients, referred to as 
the Lame and rigidity moduli, respectively. The parame- 
ter 6 G [0, 1] is the effective stress coefficient (also known 
as the Biot- Willis parameter), which describes how the 
lithostatic pressure changes with the ffuid pressure pf. 
A pressure pf > corresponds to a ffuid pressure that 
is larger than the atmospheric pressure. We use the fol- 
lowing sign conventions: > and an > in tension 
and eB < and crn < in compression. Substituting the 
stress tensor in Eq. Q into Eq. ([2| and expressing the 
strain tensor in terms of the displacement ffeld produces 
the equation of the poroelastic system. In this paper, we 
assume that the relaxation time of the elastic system is 
small compared with the time scale of the pressure evo- 
lution. We can therefore assume that for a given ffuid- 
pressure proffle, the elastic medium is always very close 
to the equilibrium state. Thus, in what follows, we will 
be concerned only with the stationary version of Eq.([2|. 

The above formalism models an isotropic system that 
does not contain any fractures. In the numerical imple- 
mentation of the problem, a simple way to model a frac- 
ture is to set the elasticity parameters equal to zero (i.e., 
G = A = 0) for the discrete volume elements containing 
a fracture (Fig.[2|. As shown by Wangen^^ this method 
of modeling the fracture correctly produces the form of 
the stress field close to the fracture tip. A more correct 
description of the stress field close to the fracture tip re- 
quires a time-dependent grid with an extra fine mesh size 
near the fracture tip. Such a detailed description is be- 
yond the scope of the present model, in which the aim 
is to provide a qualitative description of the fracturing 
process. 



B. Pressure Equation 

The equation for the fiuid-pressure is derived from the 
following continuity equation based on fiuid-mass conser- 
vatioiPl 



Here, ("(r, t) is the fiuid content added to the bulk volume, 
and QiyY^t) is a source (sink) term that arises from the 
injection (extraction) of the fiuid. The quantity vj is the 
Darcy velocity, which can be expressed in terms of the 
fiuid pressure (p/), the permeability (k), the viscosity 
(/i/), and the gravitational force (p/g, where p/ is the 
mass density of the fiuid): vj = — (k//ij)(Vp/ — P/g). 
In the most general (anisotropic) case, the permeability 
k = [kij] is a tensor. In particular, the permeability is 
strongly anisotropic in the parts of the solid that contain 
fractures. 

Let us first consider a homogeneous poroelastic 
medium that does not contain any fracture zones. In 
this case, there are two independent mechanisms that 
contribute to the fiuid content an increase in the fiuid 
pressure without an increase in the bulk strain (i.e., the 
fiuid molecules are more densely packed) , and an increase 
in the pore volume caused by an increase in the bulk vol- 
ume (while Pf remains constant). This implies that the 
fiuid content can be written a^ 



C = Spf ^bcB, 



(7) 



where is the volumetric change of the bulk volume. 
The parameter S is the constrained specific storage fac- 
tor. Thus, the time variation of ( has two contributions: 
one term that is proportional to dpf/dt and a second 
term that is proportional to dcB/dt. The characteristic 
time scale (te) of the last term strongly depends on the 
rock type, and the time scale of the first term (tp) is pri- 
marly governed by the compressibility of the fiuid. To 
compare these two time scales, let us consider an elastic 
medium that is not confined by any external forces so 
that it is free to expand/contract in response to a pres- 
sure perturbation. The response of the bulk volume is 
then given hy Scb = H~^dpf, where is the special 
poroelastic expansion coefiicient ^. With this expression 
for the volumetric response, we obtain the following ratio 
between the two time scales: 



b 



(8) 



For hard rocks with a low porosity, i.e., ^ < 0.1 and 
b < 0.4, this time ratio is small, t^/tp < 0.1^. Thus, 
for hard-rock systems, the second term in Eq. ([7| is neg- 
ligible compared with the first term, and we therefore 
disregard it in the following equations. This corresponds 
to assuming that the porous medium is incompressible. 
We also assume that the fiuid is slightly compressible. 
In this approximation, the constrained storage factor is 
5* = 0/3/, where /3/ = pj^ {dpf /dpf) is the compress- 
ibility factor under isothermal conditions. Eq. (|6| then 
becomes 



dt 



V-—{Vpf-pfg) + Qir,t). (9) 



dt 



-V-v/ + Q(r,t). 



(6) 



Next, we treat the equation describing the pressure of 
the fracture zones. In zones with fractures, a pressure 
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gradient opens the fractures and consequently modifies 
the local permeability and porosity. The time variation 
of this process is of the same order as the pressure vari- 
ation, and is a non-negligible contribution to the pres- 
sure evolution. As a consequence, the pressure equa- 
tion contains terms that couple the pressure dynamics to 
the fracture dynamics. In addition, the permeability be- 
comes anisotropic with a much higher permeability along 
the fractures than across the fractures. As described by 
Wangen^, one method to include this coupling in a nu- 
merical implementation of the problem is to assign to 
each fractured volume element (of the discretized sys- 
tem) an effective porosity that depends on the opening 
of the fracture (Fig. [2|. The effective porosity of the 
volume element i is: 



0+(l 



frL(^)' 



where 



(i) 



(10) 



(11) 



Here, ^rac(^) volume of the fracture inside ele- 

ment i, and Vq^^ is the volume of element i. The fracture 
volume of each element is found by integrating the dis- 
placement field over the fracture surface (Fig. [2| : 



denote the continuum limit of the effective porosity. Be- 
cause of the time variation of the effective porosity, the 
continuity equation for the fluid mass now becomes: 



^eff 



1 dpf 
Pf 9t 



/>eff 



dt 



- V • vf 



(12) 



which yields the following the pressure equation: 



^^'^^^-§1 + ^ = V-^(Vp/-p;g) + Q.(13) 

The fractures also modify the permeability, which be- 
comes an anisotropic second-rank tensor. In our numer- 
ical implementation of the problem, we assign to each 
volume element a permeability tensor that depends on 
the fracture direction and the fracture opening (In the 
numerical model, we allow the fracture to propagate only 
in the x and the y directions for simplicity). As an illus- 
tration, let us consider the situation shown in Fig.[2j This 
2D system has an open fracture along the y-axis, and the 
permeability is therefore larger along the y-direction than 
along the x-axis. We model this effect by introducing the 
following tensor: 



k(i) 



. k 





ii) 







(14) 



(i) 
frac 



u{t) ■ dS. 



The effective porosity becomes equal to one if the frac- 
ture fills the entire volume element, and it reduces to the 
porosity of the rock if the fracture is closed or if element 
i does not contain a fracture. 



Volume Element / 



U^+l 




Fracture volume: 

■ dS ^ |ui+i 



/ 



FIG. 2: The figure shows a system containing a fracture (the 
blue area). The fracture volume is calculated numerically 
from the displacement field in the neighboring volume ele- 
ments of the fracture. The displacement field is illustrated 
by the red arrows. The elasticity parameters of the elements 
containing parts of the fracture are equal to zero. The volume 
of each element is V^^'' = . 



Incorporating the effects of the fractures with an ef- 
fective porosity is a suitable approximation as long as 
the typical length scale of the pressure variations is large 
compared with the opening of the fracture. Let (/)eff(r, t) 



where k 



/Crock is the permeability across the frac- 



ture, which is equal to the permeability of the rock, and 



Si) 



^11 - /^rock + ( ^fluid - ^rock ) (/>frL permeability along 

the fracture. The quantity /cfluid represents the perme- 
ability inside the fracture. In the parallel plate model 
for a single fracture, the fracture permeability is given 
by the cubic law /cfluid = where w is the fracture 

aperture. This yields a very large permeability that niay 
cause numerical problems. As mentioned by Wangerv^^ 
a practical solution to this problem is to choose a perme- 
ability that is large enough to enact a minimal pressure 
drop along the fracture, but is small enough to avoid nu- 
merical instabilities. Eq. (14) becomes k = /Crockl (I is 



the identity matrix) for elements that are not fractured or 
that contain a closed fracture. The permeability tensor, 
as defined in Eq. (HS, is position-dependent and there- 



fore gives rise to terms that are proportional to V/c^j in 
Eq. lol). 



C. Fracture Criteria 

Porous materials fail under the action of a large fluid- 
pressure gradient, and the failure can be induced by both 
shear and tensile forces. Several phenomenological failure 
criteria exist^^^^. To capture both tensile and shear 
failure, we adopt two distinct criteria. 

We model the tensile failure by locally defining (i.e., for 

each volume element i) a critical tensile stress, cjc > 0. 
If one of the eigenvalues of the Terzaghi effective stress 
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tensor (evaluated at element i) exceeds the critical stress, 
the volume element i fails. The compressive strength of 
the material can be included in a similar manner. How- 
ever, during hydraulic fracturing, the tensile failure is the 
dominant fracturing mechanism, and we therefore disre- 
gard compressive failure in our numerical simulation. 

The shear failure is modeled in 2D by the Mohr- 
Coulomb criteriorP^ 



.(0 



)/2| > 



c*^*^ cos^ 



Here, is the internal friction angle, which depends on 
the density, surface, and shape of the particles constitut- 
ing the solid, and the cohesion, c, describes the minimal 
shear force that is required for fracturing when no normal 

(i) 

stresses are present. Furthermore, aa are the eigenval- 
ues of the Terzaghi effective stress tensor evaluated for 
volume element i. In a 3D system, the Drucker-Prager 
criterion is used^: 



> 



(1 



a] 



- aa, 



where r^*^ is: 

The parameter a > 1 is a dimensionless material param- 
eter. 

As discussed by Wangen^, the critical stresses (i.e., 

the parameters ai^^ and c^*^ in our model) depend on 
the grid size because the stress field is singular at the 
fracture tip. Wangen'^^ solves this problem by scaling 
the critical stress with the square root of the grid size. 
An alternative procedure, is to use experimental data 
to fit the values for ctc*^ and c^'^\ such that fracturing 
occurs for fiuid pressures in the experimental range. In 
this paper, we apply the latter approach by tuning the 
mean values of ai^^ and c^*^ so that shear and tensile 
failure appear for bore-hole pressures that are typical for 
the present problem. 

The strength of a material usually follows a distribu- 
tion that is well-described by the Weibull distribution 
(see Fig{l]i)^: 



rn I a - ath 
^0 V ^0 



exp 



Here, m is the Weibull modulus, ctq is the mean value of 
the critical stress, and ath is the threshold stress below 
which no failure will occur (usually, it is set as ath = 
0). We incorporate the heterogeneity of the rock into 
our model by distributing the local strengths ctc*^ > 
and c*^*^ > according to the Weibull distribution. The 
Weibull modulus is a measure of the degree of disorder in 
the system. A system with a small Weibull modulus has 
a higher degree of disorder than a system with a larger 
modulus. This method of modeling the heterogeneity 
is analogous to what is done in microscopic fracturing 
models, such as the fuse - and the beam -models. 



III. NUMERICAL SOLUTION APPROACH 

The model presented in the previous section is an effec- 
tive continuum model of a fractured poroelastic system. 
Our primary aim is not to provide a detailed description 
of the stress field close the fracture tip or the fiuid fiow 
inside the fractures. Instead, we seek a simplified descrip- 
tion of the coupled fluid-rock system that captures the 
primary elements of the fracturing process to obtain a 
qualitative understanding of hydrofracking, for example, 
to sort out the mechanisms that dominate the process or 
to map out what types of fracturing regimes are produced 
by different material parameters. 

To solve the system of equations presented in the previ- 
ous section, we apply a sequential solution strategy using 
standard numerical discretization methods. 



Discretization Methods 



The pressure equation, Eq. (13), is solved with a finite 



difference formulation of the complete pressure equation 
that captures the pressure dynamics in the fracture zones 
and in the homogeneous regions of the system. To iso- 
late the effect of interactions between the two injection 
points, we disregard the gravitational force in the Darcy 
velocity, and for simplicity, we consider a 2D rectangular 
region that is discretized with a lattice constant, a, as 
illustrated in Fig. [l]^. An explicit scheme is implemented 
using a forward Euler discretization in time for the tem- 
poral pressure derivative. The time derivative, d(j)ef[/dt, 
is calculated from the last two time steps. 

The ordinary differential equations arising from the 
discretized pressure equation are solved using the Cash- 
Karp embedded Runge-Kutta methocf^, which is an 
adaptive algorithm that regulates the time step using an 
error estimate. 

The boundary condition for the pressure is pj = 0. 
The Galerkin method with bilinear trial functions is used 
to solve the stationary stress equation, Eq.([2| (see Lang- 
tangen^'^ for more details). We use traction- free bound- 
aries as boundary conditions in the present problem, i.e., 
jQ^aijUjdS = for all where dfl is the boundary of 
the system, and n is its surface normal. 

Our simplified method of modeling the fractures re- 
sults in a grid-dependent fracture aperture and effective 
porosity. In practice, we fix the grid size by claiming that 
the opening of the fracture is of the order 1 mm when 
the pressure inside the fracture is of the order 1 MPa. 



B. The Fracture Event 

The present model assumes that the reservoir consists 
of hard rocks with a low porosity, so that the ratio 
in Eq. ([8| is small. Physically, this means that the 
relaxation time of the rock system is much less than 
the time variations of the pressure distribution. We 
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can therefore assume that the fracture event happens 
instantaneously. This leads to a sudden fluid-pressure 
drop inside the fracture, while the pressure outside the 
fracture is not affected. The new fluid pressure inside 
the fracture is determined using the Newton iteration, 
based on the assumption that the mass of fluid inside 
the fracture is preserved during the fracture event. The 
iterative scheme is: 



1. Adjust pf ^pf^. 

2. Solve stationary Eq. Q with the new pressure pro- 
file. 



3. Calculate the new fracture volume, V^^^^. 
iteration process. If not, return to step 1. 



Here, T) is the equation of state of the fluid, which 
expresses the fluid-mass density as a function of the tem- 
perature and the pressure. 



C. The Algorithm 

We end this theory section with a brief pseudo code 
of the numerical solution strategy: 



1. Solve the pressure equation, Eq. (13), for the cur- 
rent time step. 

2. Solve the stationary stress equation, Eq. ([2|, with 
the new pressure profile. 

3. Check the fracture criteria for each volume element. 

4. If no fracturing occurs: 

(a) Calculate the new fracture volumes. Update 

(b) Calculate a new time step from the error esti- 
mate, and return to step 1. 

5. If fracturing occurs: 

(a) Set G = A = for the fractured volume ele- 
ments. 

(b) Find the new fluid pressure and the volume 
of e ach fra cture using the iterative scheme in 



Sec. 



IIIB 



Update d^^^^/dt, and k^^). 

(c) Calculate a new time step from the error esti- 
mate, and return to step 1. Use the updated 
pressure profile as the initial pressure for the 
next time step. 



IV. MODEL SPECIFICATIONS 

Next, we apply this fracturing model to the study of 
the geomechanical stress interaction between two injec- 
tion points during the hydrofracking of a rock system 
with low permeability. 

The elasticity parameters of the system are G = 
13.7 GPa and A = 21.7 GPa. The mass-density of the 
fluid-rock system is p = 2620 kg/m^. The effective 
stress coefficient is b = 0.38, and the fluid compress- 
ibility is f3f = 5.9 X 10"^° Pa-^ The viscosity of the 
fluid, the porosity and the permeability of the rock are 
/i/ = 200 X 10-^ Pa s, ^ = 0.08 and A:rock = 0.026 mD, 
respectively. The viscosity and the compressibility of the 
fluid correspond to that of pressurized water (in the liq- 
uid phase) at a temperature of approximately 413 KP. 
These poroelastic parameters are collected from the tech- 
nical data of the Los Humeros geothermal field^. With 
these material parameters, the time ratio in Eq. (|8| is 
te/tp ~ 0.1. The permeability, /cfiuid, inside the fractures 
is four orders of magnitude larger than the rock perme- 
ability. 

We represent the reservoir as a 2D discretized system. 
To solve the pressure and stress equations, i.e., Eq. ([l3| 
and Eq. ([2|, we use a quadratic grid, as illustrated in 
Fig. [l]^, with a grid size of a = 5 m. The dimensions of 
the system range from = Ly = 150—220 m, depending 
on the separation of the injection points. 

The simulation starts with two injections points sep- 
arated by a distance, L. Initially, the system contains 
no fractures. The elasticity parameters of the injection 
elements are equal to zero. The mean values for the co- 
hesion, c, and the critical tensile stress, <Jc, are 1 MPa, 
and the internal friction angle is 40 degrees. These values 
result in tensile and shear failure for bore- hole pressures, 
typically of the order of 1 — 8 MPa. We consider systems 
with a Weibuh modulus of m G {5, 15, 30} (Fig.[l]i). The 
injection points are placed along the ?/-axis, as illustrated 
in Fig. [l] 

In the present paper, our main aim is to investigate 
how the stress interaction between two injection points 
in combination with disorder influences the hydrofrack- 
ing process. To isolate this effect, we have disregarded 
the gravitational force in our numerical implementation. 
However, in the next section, we provide a discussion of 
the consequences of the gravitational force and its inter- 
play with the stress interaction. 



V. RESULTS AND DISCUSSION 

Eq. ([5| shows that geomechanical stresses arise from 
spatial variations of the displacement field. Large spatial 
variations result in large strain and stress fields. During 
hydrofracking, the stress field is largest close to the in- 
jection point and relaxes towards zero farther away from 
the point (in the absence of gravity). Far from the injec- 
tion point, the displacement field is equal to zero. The 
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a ^ J/m3 




FIG. 3: (a) Elastic energy density of a system containing two 
injection points separated by a distance of 15 m. The bore- 
hole pressure is 5 MPa, and the injection rate is 10~^ m^/s. 
(b) Elastic energy density of a system containing two injection 
points separated by a distance of 65 m. The bore-hole pres- 
sure is 5 MPa, and the injection rate is 10~^ m^/s. Note that 
the rotational symmetry of the system is broken because the 
injection element is quadratic, i.e. the rotational symmetry 
around the injection point is reduced to quadratic symme- 
try. This is why the strain energy in (b) is not isotropically 
distributed around each injection point. 



length scale over which the displacement field relaxes to- 
ward a constant vector field is referred to as the relax- 
ation length. If two injection points are separated by a 
distance less than twice the relaxation length, the dis- 
placement field around each point is influenced by the 
presence of the neighboring point. This leads to an effec- 
tive stress interaction between the two points. 

A consequence of the geomechanical stress interaction 
is that the strain energy increases in the area between 
the two points. Fig. |3^-b shows the elastic energy den- 
sity, crijCij, stored in the elements close to the two in- 
jection points when the bore- hole pressure is 5 MPa. In 
Fig. |3^, the injection points are separated by 15 m. In 
this case, there is a significant stress interaction between 



the two points, which causes the elastic energy density 
to be largest in the region between the points. The rock 
here is under particularly strong tensile stress. In Fig.[3]3, 
the points are separated by 65 m. For such a large sepa- 
ration, the interaction between the points becomes negli- 
gible, and the elastic energy density is equally distributed 
around each of the two injection points. 

In homogeneous materials under tensile stress, a frac- 
ture propagates along the direction that causes the 
largest strain energy release rate. During hydrofrack- 
ing, we therefore expect the stress interaction between 
two injection points to mediate an effective force on the 
fractures that are created close to one of the points. This 
effective fracture force is expected to drive the fractures 
toward the neighboring point because this leads to the 
highest release rate of elastic energy. In heterogeneous 
materials, this effective force is accompanied by disorder 
effects. Whether the fracturing process is disorder-driven 
or effective force-driven, i.e., as the dominant fracturing 
mechanism, depends on the degree of disorder. 

Let us, at this point, briefly discuss the effect of the 
gravitational force. The gravitational force yields a large 
compressive stress of the order —25 MPa per distance of 
1000 m beneath the surface, which is approximately one 
order of magnitude larger than the tensile stresses pro- 
duced by the injection of fluid. The effect of the grav- 
itational force is to align the fractures along g. Thus, 
there are two distinct forces acting on the fractures: one 
caused by the gravitational field that drives the fractures 
in the vertical direction, and one arising from the ge- 
omechanical stress interaction that drives the fractures 
toward neighboring injection points. The effect of the 
gravitational field is well-known. In contrast, the effect 
of the geomechanical stress interaction is new and is the 
focus of the present paper. 

To investigate how the stress interaction influences the 
fracturing process of heterogeneous systems, we calculate 
the ensemble- averaged propagation direction, (n). The 
vector, n = {ux n^), is a unit vector that denotes the 
initial propagation direction for a fracture created from 
one of the two injection points. By definition, n = (0 1) 
points toward the neighboring point. The value of n is a 
function of the micro-state of the system. In other words, 
its value depends on the distribution of the local mate- 
rial strengths and the heterogeneity of the system. To 
map the disorder dependency of the fracturing process, 
we ensemble- average, n, by averaging over several micro- 
states. A non-vanishing (n) implies that a fracture has 
a larger probability of propagating along this direction 
than in the other directions. In this case, the fracturing 
process is strongly influenced by the stress interaction 
and the location of the neighboring injection point. In 
contrast, a vanishing (n) implies an isotropic distribution 
of n, and the governing fracturing mechanism is the dis- 
order of the system. Fig. Ilk shows (riy) as a function of 
the separation length, L, of the two injection points. The 
injection points are placed along the y-axis. (rix) is not 
shown because it is equal to zero for the systems we con- 
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Injection rate [m 7s] x 10 

FIG. 4: (a) The ensemble average, (n^), as a function 
of the separation of two injection points. {riy) is stud- 

The injection rate is 



m. 



^ m^/s. (b) The ensemble average, (n^), for the injec- 
tion rates {2 • 10~^ 2 • 10~^, 2 • 10~^} m^/s. The separation 



ied for different Weibull moduli, 
2x 10" 

2-10"^, 2-10" 
between the two injection points is 15 m, and the Weibull 
modulus is 30. In both plots, the lines are guides to the eye, 
and the ensemble average is obtained by averaging over 100 
statistical realizations of the system. 



critical correlation length, the disorder effects dominate 
the fracturing process. 

Fig.jljD shows {uy) as a function of the injection rate for 
a system in which the Weibull modulus is 30, and the two 
injection points are separated by 15 m. A decreasing in- 
jection rate results in a larger {riy) value. This is because 
of a stronger stress interaction between the two points. 
The stronger interaction arises because of a slower pres- 
sure build-up at the injection points, which leads to a 
longer time interval before the fracturing appears com- 
pared with a system with a higher injection rate. The 
pressure therefore diffuses a larger distance into the rock 
medium before the fracturing appears, which results in 
a larger relaxation length of the displacement field. The 
larger relaxation length yields a larger stress interaction 
between the points. 
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sider. As expected, {uy) decreases as the separation be- 
comes larger. This is because of a weaker interaction be- 
tween the two injection points, which results in a weaker 
effective fracture force and a reduced probability for the 
fracture to propagate toward the neighboring point. How 
quickly {uy) decays depends primarily on the degree of 
disorder. For a Weibull modulus of 30, there exists an 
effective attraction toward the neighboring point for a 
separation less than 45 m. We refer to this separation 
as the critical correlation length. If the Weibull modu- 
lus is 5, which corresponds to a highly disordered system, 
the critical correlation length decreases to approximately 
30 m. For separations less than the critical correlation 
length, the fracturing process around each injection point 
is strongly correlated with the position of the neighbor- 
ing point, and the fracturing process is governed by the 
effective fracture force. For separations larger than the 



FIG. 5: (a) The histogram shows how many times the frac- 
turing process creates a connecting fracture between the two 
injection elements. For each length, L, and Weibull modulus, 
m, 10 fracturing events are simulated, (b) The fluid pressure 
along a connecting fracture between two injection points in a 
system with a Weibull modulus of 5. (c) The fluid pressure 
along a connecting fracture between two injection points in 
a system with a Weibull modulus of 30. In (b) and (c), the 
injection points are separated by 25 m. 

The effective fracture force enhances the chance to cre- 
ate a connecting fracture network between two injection 
points. Fig. |5^ shows how many times a number of frac- 
turing processes creates a connecting fracture between 
the two injection points as a function of the separation for 
systems with different degrees of disorder. Two examples 
of a connecting fracture are shown in Fig. and Fig. 
for systems with Weibull moduli of 5 and 30, respectively. 
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In the typical time evolution of such a fracture, the frac- 
ture first connects the two points before a new fracture 
arm is created from each point. In addition, for a system 
with a high Weibull modulus (Fig. [sj)), the fracture is 
most likely to show the same character as a fracture in a 
homogeneous system, i.e., a straight line fracture, while 
the fractures in highly disordered systems (Fig. ^) prop- 
agate more randomly. The ability to create a connecting 
fracture network is crucial for creating a geothermal sys- 
tem, such as an Enhanced Geothermal System (EGS), as 
well as in the exploitation of unconventional hydrocarbon 
resources. In the case of hydrofracking, the effect of the 
stress interaction is important. During hydrofracking, 
the perforation zones ideally are completely uncorrelated 
to create long and deep cracks that connect a large area 
of the reservoir to the well. If the zones are spaced by a 
distance of less than the critical correlation length, the 
stress interaction leads to a reduced hydraulic stimula- 
tion of the system. We therefore believe that our results 
can be used in the optimization of the hydrofracking pro- 
cess. 



VI. CONCLUSIONS 

In this paper, we developed a model of hydrofracking 
and applied the formalism to the study of how geome- 
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chanical stress interactions between two injection points 
influence the fracturing process. We found that when 
the separations between the two injection points is less 
than a critical correlation length, the fracturing process 
around each injection point is strongly correlated with 
the position of the neighboring point. The magnitude 
of the critical correlation length depends on the degree 
of heterogeneity of the rock. For weakly disordered sys- 
tems, the correlation length can be as large as 45 m, and 
for highly disordered rock systems, it reduces to approx- 
imately 20 m. In the strongly correlated regime, there 
exists an effective fracture force that drives the fractures 
toward the neighboring injection point. An important 
observation in this work is that the fracture force re- 
duces the effectiveness of the hydraulic stimulation if the 
injection points are separated by a distance less than the 
critical correlation length. 
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